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Diffusion processes and coalescent trees 

Robert C. Griffiths^ and Dario Spanq^ 



,— j- ' We dedicate this paper to Sir John Kingman on his 70th Birthday. 

In modern mathematical population genetics the ancestral history of 
a population of genes back in time is described by John Kingman's 
coalescent tree. Classical and modern approaches model gene frequencies 
Ph ' by diffusion processes. This paper, which is partly a review, discusses 

C"| . how coalescent processes are dual to diffusion processes in an analytic 

"j^ ' and probabilistic sense. 

Bochner (1954) and Gasper (1972) were interested in characterizations 
of processes with Beta stationary distributions and Jacobi polynomial 
eigenfunctions. We discuss the connection with Wright-Fisher diffusions 
J> ■ and the characterization of these processes. Subordinated Wright-Fisher 

diffusions are of this type. An Inverse Gaussian subordinator is inter- 
esting and important in subordinated Wright-Fisher diffusions and is 
■^j- ■ related to the Jacobi Poisson Kernel in orthogonal polynomial theory. 
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A related time-subordinated forest of non-mutant edges in the Kingman 



' coalescent is novel. 
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R. C. Griffiths and D. Spano 
1 Introduction 



The Wright -Fisher diffusion process {X(t),t > 0} models the relative 
frequency of type a genes in a population with two types of genes a 
and A. Genes are subject to random drift and mutation over time. The 
generator of the process is 

c = \<i-*)£ + \(-°* + W-*))b (i-i) 

where the mutation rate A — > a is \a and the rate a — > A is \fi. If a 
and /3 are zero then zero and one are absorbing states where either A or 
a becomes fixed in the population. If a, (3 > then {X(t),t > 0} is a 
reversible process with a Beta stationary density 

UAv) = B(a, /3rV _1 (l - yf~\ < y < 1. (1.2) 
The transition density has an eigenfunction expansion 

f(x,y;t) = / a>j9 (y) { 1 + f] p e „ (i)P>« (x)P>« (») ) , (1.3) 

^ n=l - 1 

where # = a + /3, 

p£(t)=exp{-in(n + 0-l)*}, (1.4) 

and {P^ Q,,9 ' ) (y), n G are orthonormal Jacobi polynomials on the 

Beta (a, /?) distribution, scaled so that 



E 



P^\Y)P^\Y) 



m, n G Z + 



under the stationary distribution (|1 .2|) . The Wright-Fisher diffusion is 
also known as the Jacobi diffusion because of the eigenfunction expansion 
(jl.3p . The classical Jacobi polynomials, orthogonal on 

(l-x) a {l + xf, -l<x< 1, 

can be expressed as 

P K/3) (a; ) = ( a + 1 \n) 2Fli _ njn + a + ^ + i ;a + 1; (i _ x)/2), (1.5) 

where 2 Pi is a hypergeometric function. The relationship between the 
two sets of polynomials is that 

P^\x) = c n P^- 1 ' a - 1 \2x-l), 
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where 



Define 



(2n + a + (3 - l)(a + /3) (n _ 1) n! 

«(n)/3(n) 



the forward generator of the process. The Jacobi polynomials are eigen- 
functions satisfying, for n € Z+, 

CP^\x) = -\n{n + e-l)P^\x); 

£f a ,p{x)pW{x) - -±n(n + 9 l)f ai p(x)P^(x). (1.7) 

The well known fact that the Jacobi polynomials {Ph a '^(x)} satisfy 
(jl.7l) implies that they are eigenfunctions with corresponding eigenvalues 

{/*(*)}■ 

In modern mathematical population genetics the ancestral history of a 
population back in time is described by John Kingman's elegant coales- 
cent process jl9^ . The connection between the coalescent and Fleming- 
Viot diffusion processes is made explicit by Donnelly and Kurtz in [7J| , 
Q by their look-down process. An approach by Ethier and Griffiths fioj j 
uses duality to show that a 'non- mutant lines of descent' process which 
considers a forest of trees back in time to their first mutations is dual 
to the Fleming-Viot infinitely-many-alleles diffusion process. The two- 
allele process {X(t),t > 0} is recovered from the Fleming-Viot process 
by a 2-colouring of alleles in the infinitely-many-alleles model. If there 
is no mutation then the dual process is the same as the Kingman co- 
alescent process with an entrance boundary at infinity. The dual process 
approach leads to a transition density expansion in terms of the trans- 
ition functions of the process which counts the number of non-mutant 
lineages back in time. It is interesting to make a connection between 
the eigenfunction expansion (ll.3[) and dual process expansion of the 
transition densities of {X(t),t > 0}. Bochner Q and Gasper [li| find 
characterizations of processes which have Beta stationary distributions 
and Jacobi polynomial eigenfunctions. Subordinated Jacobi processes 
{X(Z(t)j ,t > 0}, where {Z(t),t > 0} is a Levy process, fit into this 
class, because subordination does not change the eigenvectors or the 
stationary distribution of the process. The subordinated processes are 
jump diffusions. A particular class of importance is when {Z(t),t > 0} 
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is an Inverse Gaussian process. Griffiths [l8[ obtains characterizations 
of processes with stationary distributions in the Meixner class, as well 
as for Jacobi processes. The current paper is partly a review paper de- 
scribing connections between Jacobi diffusions, eigenfunction expansions 
of transition functions, coalescent trees, and Bochner characterizations. 
Novel results describe the subordinated non-mutant lines-of-descent pro- 
cess when the subordination is with an Inverse Gaussian process. 



2 A coalescent dual process 



A second form of the transition density (|1 .3[) derived in Ethier and Grif- 
fiths is 



f(x,y i t) = ^^)^2B(l-,k,x)f a+ i, 0+k - l { 3 f), (2.1) 



k=0 1=0 



where 



B(l;k,x) = ( )x fc (l-x)^ fc , fc = 0, 



is the Binomial distribution and {<Z/?(£)} are the transition functions of 
a death process with an entrance boundary of infinity, and death rates 
k(k + 9 — I)/2, k > 1. The death process represents the number of non- 
mutant ancestral lineages back in time in the coalescent process with 
mutation. The number of lineages decreases from k to k — 1 from co- 
alescence at rate Q) or mutation at rate k0/2. If there is no mutation, 
{qu(t),t > 0} are transition functions of the number of edges in a King- 
man coalescent tree. There is an explicit expression for the transition 
functions beginning with the entrance boundary of infinity 16, 2lL 
of 



9fc(*) = E'°iW(- 1 > 



,._ fc (2j+g-l)(fc + %_ 1) 
k\{j-k)\ 



(2.2) 



recalling that p s n {t) is defined by (11.41) . A complex- variable representation 
of (I2.2[) is found in Let {X t ,t > 0} be standard Brownian motion 
so X t is N(0, t). Denote Z t = e iXt and u> t = e^^ et , then 



q ^>- e8 T(k + 9)kl E 



(u> t Z t ) k (l-uj t Z t ) 



rz~ t (\ + uj t z t y k +° 



(2.3) 
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for k = 0, 1, .... The transition functions for the process beginning at 
n, rather than infinity, are 



j—k 



for k = 0, 1, . . . , n. An analogous complex-variable representation to 
(f273|) is 



?nfcw T{k + 6)r(n + k + 6)\k 



x 2 i ? i(-n + fc + 1, + 2fc; n + fe + 0; Z t )] (2.5) 

for k = 0, 1, . . . , n. The expansion (|2.1[) is derived from a two-dimension- 
al dual death process {L 6 (t) e Z 2 + ,t > 0} which looks back in time in 
the diffusion process {X(t),t > 0}. A derivation in this paper is from 
, which follows more general analytic derivations in [lol ] for a Fleming- 
Viot model and Q for a diffusion model with selection. Etheridge and 
Griffiths Q give a very clear probabilistic derivation in a Moran model 
with selection that provides an understanding of earlier derivations. A 
sketch of a derivation of (|2.ip from [9( is the following. Let x\ = x, 
x-2 = 1 — x and define for k £ 

/ \ ^(l fc l) fci fe 2 
9k(x) = x 1 1 x 2 a , 

a (ki)P(k 2 ) 

then 

£9k(x) = \{\k\+9- l)[k l9k - ei (x) + k 2 g k _ e2 (x) - \k\g k (x)]. (2.6) 

Here and elsewhere we use the notation \y\ = ^ - =1 Uj for a d-dimension- 
al vector y. In this particular case \k\ = k% + k 2 . To obtain a dual process 
the generator is regarded as acting on k = (ki,k 2 ), rather than x. The 
dual process is a two-dimensional death process {L (t),t > 0}, the rates 
of which are read off from the coefficients of the functions g on the 
right-hand side of (|2.6[) ; 

1 k- 

k^k-Bi at rate --^ • Ifclflfcl + 9 - I). (2.7) 

2 |fc| 

The total size, \L (t)\, is a I-dimensional death process in which 

at rate ~\k\(\k\ + - 1) 
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with transition functions denoted by {q e m i{t),t > 0}. There is hypergeo- 
metric sampling of types which do not die, so 



P(L(t) = I | L(0) = m) = q e ml (t) = qf mm (t) 



(T)( 



m 2 \ 



GDll'l ' 



(2.8) 



where q\ m \\i\{t) is defined in (|2.4|) . The dual equation obtained by re- 
garding L as acting on x or k in (|2.6p is 



E 



X(0) 



<?L(0)(*W) 



= E 



L(0) 



5L( t )(X(0)) 



(2.9) 



where expectation on the left is with respect to the distribution of X(t), 
and on the right with respect to the distribution of L(t). Partitioning 
the expectation on the right of (|2.9I) by values taken by L(t), 



E, 



x 1 (t) m ^x 2 (ty 



\m\ 
mi 

m|\ Q!(mi)^(m 3 ) 



™1/ ^(rai+m 2 ) 



Ed r 
2-1 ^2 



'(1*1) 



Km 



(2.10) 

Z|\ mi [;i] TO 2[/2] 



/i 



l TO ki'i] 



The transition distribution of X(t) now has an expansion derived from 
an inversion formula applied to (|2.10p . Letting mi, m 2 — > oo with 
mi/\m\ -> 2/1, m 2 /|m| -> y 2 gives 



/e: 



/(«, j/; t) = J2 (*) ( i x ) aM^a + «x,y9 + Ja)" 1 ^ 



l„,Ji+a-l,Ja+i3-l 
i/2 > 



which is identical to (|2.1j) . 

The two-allele Wright -Fisher diffusion is a special case of a much more 
general Fleming- Viot measure- valued diffusion process which has V(S), 
the probability measures on S, a compact metric space, as a state space. 
The mutation operator in the process is 

(Af)(x)= 9 - f (/(£)- /(z)H(^), 
J s 

where vq S V(S) and / : 5 —> HL The stationary measure is a Poisson- 
Dirichlet (Ferguson-Dirichlet) random measure 

oo 

A* = X^ l< ^' 

i=i 

where {xi} is a Poisson-Dirichlet point process, VT>(6), independent of 
{£j} which are i.i.d. v £ V(S). A description of the VT>(9) distribution 
is contained in Kingman [20j |. 
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Denote the stationary distribution of the random measure as 

IV (-)=P(/iG-)- 

Ethier and Griffiths [loj ] derive a transition function expansion for P(t, fj,, 
dv) with given initial /i € V(S) of 

p(t,Ai,.) = 8g(t)n fl , Wb (0 

+ E li(t) / x • • • x dx n ) 

n=l J S" 

x n n+ i ( n+ e)-i{„^ n ( Xli ... )a . n ) + I/o }.(-), (2.11) 

where r) n (xi, . . . , cc„) is the empirical measure of points x\, . . . , x n G S: 

r/ n (xi,...,x n ) = n~ 1 (6 Xl H \-S Xn ). 

There is the famous Kingman coalescent process tree [l9| behind the 
pretty representation (|2.1ip . The coalescent tree has an entrance bound- 
ary at infinity and a coalescence rate of (fj) while there are k edges in 
the tree. Mutations occur according to a Poisson process of rate 9/2 
along the edges of the coalescent tree. (<z^(i)} is the distribution of the 
number of non-mutant edges in the tree at time t back. The number of 
non-mutant edges is the same as the number of edges in a forest where 
coalescence occurs to non-mutant edges and trees are rooted back in 
time when mutations occur on an edge. If the time origin is at time t 
back and there are n non-mutant edges at the origin then the leaves of 
the infinite-leaf tree represent the population at t forward in time di- 
vided into relative frequencies of families of types which are either the 
n non-mutant types chosen at random from time zero, or mutant types 
chosen from vq in (0,i). The frequencies of non- mutant families, scaled 
to have a total frequency unity, have a Dirichlet distribution with unit 
index parameters, and the new mutation families, scaled to have total 
frequency unity, are distributed according to a Poisson-Dirichlet random 
measure with rate 9 and type measure Uq. The total frequency of new 
mutations has a Beta {9, n— 1) distribution. An extended description of 
the tree process is in Griffiths (l7j . 

A (i-dimensional reversible diffusion process model for gene frequencies 
which arises as a limit from the Wright-Fisher model has a backward 
generator 

* = \ EE ^ - + \ X> - (2-12) 

2—1 7 — 1 i—1 
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where 9 = \e\. In this model mutation is parent-independent from type 
i — > j at rate ^6j, i, j — 1, . . . , d. Assuming that e > 0, the stationary 
density is the Dirichlet density 



^-•••xT> ( 2 - 13 ) 



r( ei )-..r(e d ) 1 

for xi, . . . , Xd > and Xi — 1. Griffiths [15| shows that the transition 
density in the model has eigenvalues 

P|n| (t) = e -iH(N+ e - 1 )* 

repeated 

\n\ +d- 2' 
\n\ . 

times corresponding to eigenvectors {Q° (x), n g Z'fT 1 } which are multi- 
type orthonormal polynomials of total degree |n| in x. As eigenfunctions 
the polynomials satisfy 

£Q° n {x) = -\\n\{\n\ + 9 - l)Q° n {x). (2.14) 

The eigenvalues {pk{t), k £ Z + } do not depend on the dimension d. The 
transition density with X(0) = x, X(t) = y has the form 

f(x,y,t)=V(y,e)il+ p H (t)Q H (x,y)\ . (2.15) 

The kernel polynomials on the Dirichlet {Q\ n \ (x, y)} appearing in (|2.15p 
are defined as 

Q H (x,y)= Qn( x )Q°n(y) (2-16) 

{n:|n| fixed} 

for any complete orthonormal polynomial set {Qn( x )} on the Dirichlet 
distribution ([2~T31) . If d = 2, 

Q H (x,y)=P^\x)P^\y) 

where {P\n\' €2 \ x )} are orthonormal Jacobi polynomials on the Beta 
distribution on [0, 1]. In general n is just a convenient index system for 
the polynomials since the number of polynomials of total degree |n| is 
always the same as the number of solutions of ni + •• • + n c i-i = \n\, 

n\+d-2 
\n\ 
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Q\ n \ i x t y) is invariant under the choice of which orthonormal polynomial 
set is used. The individual polynomials Qn(x) are uniquely determined 
by their leading coefficients of degree \n\ and Q\ n \(x,y). A specific form 
is 

Q H (x,y) = (9 + 2\n\ - 1) ]T (-l)H-m ( g + H(N-i) ^ (2 1?) 

* — ' to n — to) 

m=0 

where 



|;|=m 

An inverse relationship is 



6»= E (2-18) 



The transition distribution f|2 . 15[) is still valid if any or all elements of 
e are zero. The constant term in the expansion then vanishes as the 
diffusion process is transient and there is not a stationary distribution. 
For example, if e = 0, 

d ( oo ^ 

f{x,y,t) = Y[yf\ E P\n\(t)Cfl n{ (x,v) , (2.20) 
3=1 l\n\>d J 



where 



Q^(x,y) = (2M - i)^(-i)N-m ^(W-i) «o 
11 ^-^ m!(|n| — to)! 



with 



{/:/>0,|/|=m} V 7 111 V** 1 >- 1 

The derivation of (|2.15[) is a very classical approach. The same process 
can be thought of as arising from an infinite-leaf coalescent tree similar to 
the description in the Fleming- Viot infinitely-many-alleles process. The 
coalescent rate while there are k edges in the tree is and mutations 
occur along edges at rate 8/2. In this model there are d types, 1,2, . . . , 
d and the probability of mutation i — > j, given a mutation, is ej/9. This 
is equivalent to a d-colouring of alleles in the Fleming- Viot infinitely- 
many-alleles model. Think backwards from time t back to time 0. Let 
y = (yi, . . . , yd) be the relative frequencies of types in the infinite number 
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of leaves at the current time t forward and x — (x\, . . . , X4) be the 
frequencies in the population at time 0. Let / be the number of non- 
mutant edges at time which have families at time t in the leaves of the 
tree. Given these I edges let U = {U\, . . . , Ui) be their relative family 
sizes in the leaves, and V — (Vi, . . . , Va) be the frequencies of families 
derived from new mutations on the tree edges in (0, t). The distribution 
of U(SV = (Ui,...,U h Vx,...,V d ) is23(«e»,(l,...,l)®e). The type 
of the I lines, and therefore their families, is chosen at random from the 
frequencies x. The distribution of the number of non-mutant lines at 
time from the population at t is qf(t). The transition density in the 
diffusion (|2.15[) is identical to the mixture distribution arising from the 
coalescent 

00 

f(x,y,t)=J24\( t ) E M(l,x)V(y,e + l), (2.23) 

|i|=0 {l:\l\ fixed} 

by considering types of non-mutant lines, and adding Dirichlet variables 
and parameters according to U non-mutant families being of type i. 
M.{1, x) is the multinomial distribution describing the choice of the initial 
line types from the population at time 0. The expansion when d = 2 
corresponds to (|1.3[) . The argument is valid if any elements of e are zero, 
considering a generalized Dirichlet distribution T>(x, e) where if e; = 0, 
then Xi — with probability 1. 

The algebraic identity of (|2.23l) and (I2.15[) is easy to see by expressing 
Q\ n \(x, y) in terms of {£ m }, then collecting coefficients of in (|2.15[) to 
obtain (|2.23|) . Setting po(t) = 1 and Qo{x, y) = 1, the transition density 
is 

CO 

f(x,y,t)=V(y,e) p\ n \(t)Q\ n \(x,y) 

|n|=0 

00 

= E 

x V(y,e)£i(x,y) 

OO 

= E#) E M(l,x)V(y,e + l). (2.24) 

1*1=0 {l:\l\ fixed} 

The non-mutant line-of-descent process with transition probabilities 
{?n(*)} appears in all the Wright-Fisher diffusion processes mentioned 
in this section as a fundamental dual process. The process does not de- 
pend on the dimension of the diffusion, partly because the (/-dimensional 



n\=\l\ 



..,(t)(0 + 2|n|-l)(-l) 



\n\-\l\ 



(0 + 

|I|!(|n|-|Z|) 



){\n\--q_ 

! 
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process can be recovered from the measure-valued process as a spe- 
cial case by colouring new mutations into d classes with probabilities 
(ei/6, €2/0, . . . , €4/6) with 6 = Y^j=i € j- ^ * s a l so interesting to see the 
derivation of the d-dimensional transition density expansion as a mix- 
ture in terms of {q^t)} via the orthogonal-function expansion of the 
transition density in (|2.24p . 



3 Processes with beta stationary distributions and 
Jacobi polynomial eigenfunctions 

In this section we consider 1-dimensional processes which have Beta sta- 
tionary distributions and Jacobi polynomial eigenfunctions, and their 
connection with Wright-Fisher diffusion processes. We begin by consid- 
ering Bochner Q and Gasper's [l3| characterization of bivariate Beta 
distributions. 

A class of bivariate distributions with Beta marginals and Jacobi poly- 
nomial eigenfunctions has the form 

f(x,y) = f a0 {x)f a p{y)ll + f^p n P^\x)P^\y)V (3.1) 
^ n=l ' 

where {p„, n S Z + } is called a correlation sequence. The transition dens- 
ity (|1.3p in the Jacobi diffusion has the form of the conditional density 
of Y given X = x in (|3.1[) with p n = p e n {t). Bochner Q and Gasper [l3[ 
worked on characterizations of sequences {p n } such that the expansion 
(|3.1I) is positive, and thus a probability distribution. It is convenient to 
normalize the Jacobi polynomials by taking 

so that Ri^' 13 ^!) — 1; denote 

a {n) n\ 

and write 

f{x,y)=f a0 {x)f a ^{l + f^p n KRt P \x)R^\y)\. (3.2) 
^ n=l ' 

Bochner Q defined a bounded sequence {c„} to be positive definite with 
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respect to the Jacobi polynomials il 

^a n h n R^(x)>0, I 

a n \h n < oo 



n>0 



n>0 



implies that 



^a nCn h n Ri a <P\x) >0. 



n>0 



Then is a correlation sequence i/ and on/?/ if it is a positive definite 
sequence. The only if proof follows from 



J2 a nPn Rt P) (x)=E a n h n Ri a ^ (Y) 



n>0 



- ri>0 



X = x 



>0, 



where (X, Y) has the distribution (|3.2j) . The if proof follows at least 
heuristically by noting that 



n>0 



where S(-) has a unit point mass at zero, so if {p n } is a positive definite 
sequence then 

Y / PnKR^\x)Ri a ^{y)>Q 



n>0 



141. 



and (|3.2p is non-negative. A careful proof is given in 
Under the conditions that 



a < j3 and either 1/2 < a or a + j3 > 2, (3.3) 

it is shown in [l3| that a sequence p n is positive definite if and only if 

p n = E[R^(Z)] (3.4) 

for some random variable Z in [0, 1]. If the conditions (|3.3[) do not hold 
then there exist x, y, z such that K(x, y, z) < 0. The sufficiency rests on 
showing that under the conditions (13.31) for x, y, z G [0, 1], 



K(x,y,z) = J2h n Rt^(x)Rt^(y)R^(z) > 0. 



(3.5) 



The sufficiency of (|3.4[) is then clear by mixing over a distribution for Z 
in (|3.5[) to get positivity. The necessity follows by setting x = 1 in 

and recalling that i?£ Q '' 3 ' ) (1) = 1, so that Z is distributed as Y conditional 
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on X = 1. This implies that extreme correlation sequences in exchange- 
able bivariate Beta distributions with Jacobi polynomial eigenfunctions 
are the scaled Jacobi polynomials {R%*'^\z), z <E [0,1]}. Bochner 0] 
was the original author to consider such problems for the ultraspherical 
polynomials, essentially orthogonal polynomials on Beta distributions 
with equal parameters. 

A characterization of reversible Markov processes with stationary Beta 
distribution and Jacobi polynomial eigenfunctions, from [13] . under 
|, is that they have transition functions of the form 



f(x,y;t) = f a0 (y)i 1 + £ c n {t)h n R^ {x)R^ {y) 

^ n=l 

with c n (t) — exp{— d n t}, where 



d* 



an(n + a + j3 - 1) + 



1-z 



v(dz) 



(3.6) 



(3.7) 



g > 0, and v is a finite measure on [0, 1). If v{ ) = 0, a null measure, 
then f(x,y;t) is the transition function of a Jacobi diffusion. 

Eigenvalues of a general reversible time-homogeneous Markov process 
with countable spectrum must satisfy Bochner's consistency conditions: 

(i) {c n (t)} is a correlation sequence for each t > 0, 

(ii) c n (t) is continuous in t > 0, 

(iii) c n (0) = c (i) = 1, and 

(iv) c n (i + s) = c„(i)c„(s) for t, s > 0. 

If there is a spectrum {c„(i)} with corresponding eigenfunctions {£ n } 
then 

c n (t + s)^ n (X(0))=E^ n (X(t + s)) | X(0) 

= E[E[e„(X(f + s)) | X(s)] X(0) 

= c n {t)E^ n (X(s)) \ X(0) 

= c n (t)c n (s)Z n (X(0)), 

showing (iv). If a stationary distribution exists and X(0) has this distri- 
bution then the eigenfunctions can be scaled to be orthonormal on this 
distribution and the eigenfunction property is then 



£ m (X(t))£ n (X(0)) =c n (t)6 
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{X(t),t > 0} is a Markov process such that the transition distribution 
of Y = X(t) given X(0) = x is 

f(x, y; t) = f(y)!l + ]T c n (t)U^)Uy)\, (3-8) 
^ n=l ' 

where /(y) is the stationary distribution. In our context are the 

orthonormal Jacobi polynomials. A Jacobi process {X(t),t > 0} with 
transition distributions p.6[) can be constructed in the following way, 
which is analogous to constructing a general Levy process from a com- 
pound Poisson process. Let {Xk,k € Z + } be a Markov chain with sta- 
tionary distribution f a p (y) and transition distribution of Y given X = x 
corresponding to (13. 1|) . with (|3.3I) holding, and {N(t),t > 0} be an in- 
dependent Poisson process of rate A. Then (Xo,Xk) has a correlation 
sequence {p 1 ^} and the transition functions of X(t) — X^m have the 
form fljTSJ), with 

d n = X f (1-Rl a ^(z))n(dz), (3.9) 
Jo 

where (i is a probability measure on [0,1]. The general form (|3.7p is 
obtained by choosing a pair (A, /ia) such that 

r 1 r 1 1 - R^'^M 
d„ = lim A / (l-i4 a,fl) («))MA(«te)= / —v{dz). (3.10) 

Equation (|3. 101) agrees with p.7[) when any atom ^({1}) is taken out of 
the integral, because 

lim \li = cn ( n + 6-1), 

z— s-1 1 — Z 

where c > is a constant. 



4 Subordinated Jacobi diffusion processes 

Let {X(t),t > 0} be a process with transition functions p.6p . and 
{Z(t),t > 0} be a non-negative Levy process with Laplace transform 

-tj H{dy)\, (4.1) 

where A > and H is a finite measure. The subordinated process 
{X(t) = X(Z(t)),t > 0} is a Markov process which belongs to the 



e- d ^G{dy) 
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same class of processes with correlation sequences 

c n {t) = E[c n (Z(t))] = exp |-t jf H{dy)\ , (4.2) 

where _ff is a finite measure. c n (t) necessarily has a representation as 
e~ dnt , where d n has the form (|3 . 10[) for some measure v. We describe 
the easiest case from which the general case can be obtained as a limit. 
Suppose 

Jo V 

and write 

ff(riy) 

G = — : — , 
Ay 

so that G is a probability measure. Let 

K(dz) = f aP (z)dz!l + f]h n R^(z) f 
Then if is a probability measure and 

A jf (l - R^\z))K(dz) = A jH(l - e^)G(dy). 

The representation (|3.10[) is now obtained by setting 

= A(l - z)K(dz). 

We now consider subordinated Jacobi diffusion processes. The subor- 
dinated process is no longer a diffusion process because {Z(t),t > 0} 
is a jump process and therefore {X(t),t > 0} has discontinuous sample 
paths. It is possible to construct processes such that (|4.2p holds with 
d n = n by showing that e~ tn is a correlation sequence and thus so is 
E[e — ^W n l. The construction follows an idea in Q. The Jacobi-Poisson 
kernel in orthogonal polynomial theory is 

oo 

l + J2r n h n Ri a ^(x)R^(y) 7 (4.3) 
n=l 

which is non-negative for all a, (3 > 0, x, y € [0, 1], and < r < 1, for 
which see pi 12. The series (I4.3p is a classical one evaluated early in 
research on Jacobi polynomials (see Q). In terms of the original Jacobi 
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polynomials, (|1.5[) 



n=0 



r(g + /3 + 2)(l -r) 

2«+/ ? + 1 r(a + l)r(/3 + 1)(1 + r) a +P+ 2 

((a + P + 2)/2) (m+n) ((a + g + 3)/2) (m+n) / ffl 2 
(a + l) (m) (/3 + l) (m) m!n! U 2 

, 2 Q +^ +1 r(n + a + l)r(n + ^ + l 



E 



fc 2 J 



(4.4) 



where 



2n + a + /3 + 1 T(n + l)T(n + a + /3 + 1) ' 



£ = cos2iy9, y = cos 20, a = sin sin 0, 6 = cosy cos 0, k — (r 1 / 2 + 
r~ 1/2 )/2. The series (|4~4l) is positive for -1 < x, y < 1, < r < 1 and 
a, j8 > -1. 

A Markov process analogy to the Jacobi-Poisson kernel is when the 
eigenvalues c n (t) — exp{— nt}. Following jfjj let X(t) = X{Z{£j), where 
{Z(t),t > 0} is a Levy process with Laplace transform 



E 



,-\z{t) 



exp 



exp 



{-t[y/2X + (6- l) 2 /4 - y/(0 - l) 2 /4] } (4.5) 



{Z(t),t > 0} is a tilted positive stable process with index ^ such that 
Z(t) has an Inverse Gaussian density 

2t 
.\0~ 

that is, 

t 



# i; 



\/2irz- 



: exp 



1 /l^- 1| \ 2 1 

' — z- t) j, z > 0. 



(4.6) 



The usual stable density is obtained when = 1 and (|4.6[) is a tilted 
density in the sense that it is proportional to exp{ — z{6 — l) 2 /8} times 
the stable density. See jl2] XIII, §11, Problem 5 for an early derivation. 
Z(t) is distributed as the first passage time 



T t = infju > 0;B(u) 



-u = t 
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The process {X(t),t > 0} is a jump diffusion process, discontinuous at 



the jumps of {Z(t),t > 0}. Jump sizes increase as 9 decreases. If 9 < 1 
then for n > 1 



so subordination does not directly produce eigenvalues e~ nt . Let f(x, 
y;t) be the transition density of X(t), then the transition density with 
eigenvalues exp{— nt}, n > is 



The subordinated process with this transition density is X(Z(t)), where 
Z(t) is a similar process to Z(t) but has an extra state infinity. Z(t) is 
killed by a jump to infinity at a rate (1—9). Another possible construction 
does not kill the process X, but restarts it in a stationary state drawn 
from the Beta distribution. It is convenient to use the notation that 
a process {Z°(t),t > 0} is {Z(t),t > 0} if 6 > 1, or {Z(t),t > 0} 
if < 9 < 1, and use the single notation {X(Z°(t)),t > 0} for the 
subordinated process. The transition density (|3.6p . where c n (t) has the 
general form 



can then be obtained by a composition of subordinators from the Jacobi 
diffusion with any a, > 0. 

There is a question as to which processes with transition densities 
(|3.6I) and eigenvalues c n (t) described by (|3.7[) are subordinated Jacobi 
diffusion processes. We briefly consider this question. Substituting 



E exp 




e-^- e )J(x,y;t)+(l-e-^)f a p(y). 




R { n' P) {y)= 2F 1 (-n,n + e-l;P;l-y) 
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in the eigenvalue expression (|3.7p . 



I 



v(dy) 



/c=i 



1-1/ 

{~n) {k) (n + 6 - l) (fe) ^ fc _! 



/3(fc) 



fc! 



a n?=o (-»(» + g - 1) + jq- + g - 1)) 



fc=i 



/3( fe ) 



where (l — y) k v{dy) = c^k- The generator corresponding to a process 
with these eigenvalues is 



n- ~o(2£+j(i+g-i)) 

/3( fc ) k\ ' 



where £ is the Jacobi diffusion process generator (jl.ip . The structure of 
the class of stochastic processes with the generator C needs to be un- 
derstood better. It includes all subordinated Jacobi diffusion processes, 
but it seems to be a bigger class. A process with generator £ is a sub- 
ordinated Jacobi diffusion process if and only if the first derivative of 

% fcT (4 ' 8) 

is a completely monotone function of A. Factorizing 

-2A + j(j +6-1) = (j + n (A))(j + r 2 (A)), 
where ri(A),r2(A) are 



fc=i 



1)/2±V2A + (0-l) 2 /4, 



(|4~8l) is equal to 
i-i- 



o 



2 F 1 (r 1 (A),r 2 (A);/3; 1 - y) - 1 (1 - y)-^(dy). 



(4.9) 



5 Subordinated coalescent process 

Subordinating the Jacobi diffusion process {X(t),t > 0} leads to subor- 
dinating the coalescent dual process, which we investigate in this section. 
A subordinated process {X(t) — X(Z(t)),t > 0} has a similar form for 
the transition density as (|2.ip . with qf(t) replaced by E(qf(Z(t)), which 



Diffusion processes and coalescent trees 19 

are transition functions of the subordinated death process A e (Z(t)). The 
subordinated process comes from subordinating the forest of non-mutant 
lineages in a coalescent tree. 

If A e (t) = A 6 (Z°(t)), with Z°[t) defined in the last section, we will 
show that the probability distribution of A e (t) , 9 > is 

{ k )(— ) (— ) ^ ^ 

for k € Z+, where z = e~*. The distribution (|5.ip is the distribution of 
the number of edges in a time-subordinated forest. Note that if < 9 < 1 
we still invoke a subordinator with a possible jump to infinity at rate 
1-9, so 

E[4(Z°(t))] = e-^E[4(Z(t))] + (1 - e-^)6 k0 , 

because q^(oo) = 8ko- Although 9 is greater than zero in (|5.1I) . it is 
interesting to consider the subordinated Kingman coalescent with no 
mutation. Then A°(t) > 1, and 

E[q° k (Z°(t))] = e-*E[gg(Z(i))] + (1 - e -*)* fcl , 

because a jump to infinity is made at rate 1, and <z£,(oo) = Ski- The 
distribution of A°(t) is then, for k > 1, 

(VXlT^TT^-^Mi-*). OU) 

The proof of (pTTj) (0 > 0) and (with 9 = 0) follows directly from 

the expansion (|2.2I) . 

« M( z- W )]-f^(-i)^ 3^ag±^ 



r(2fc 



fc!r(fc + 6>) 

x {l + E(-lF(2i + 2* + - l) (2fc + ^-V j 
•> j=i J 

-ffi^«*< 1 -'>< 1+ '>-'"** 1 

= ( * J (it-;) (ttJ (5 - 3) 

Effectively, in the expansion (I2.2[) of <?f (t), terms /»/(<) = exp{— + 
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8 — l)t} are replaced by z J = exp{— jt}. The third line of (|5.3p follows 
from the identity, with \z\ < 1 and a = 2k + 8, that 



(1 -z)(l + z)- a = 1 + X>l)'(2j + a - 1) 



"O-D -7 



3=1 J 



proved by equating coefficients of on both sides. Of course, for any 
\z\ < 1, since (|5.1[) is a probability distribution, 

2k + e-l\/ z \ k / 1 \ k + e 



fe=0 



* (— ) (1 -^ ) = L (5 - 4) 



The probability generating function of (|5.1[) is 

where p = e - */ (1 + e - *) and q = 1/(1 + e _t ). The calculation needed to 
show (15.51) comes from the identity 



V ( 2fc + * " V = ^£±^3 , (5.6) 



fc=0 

which is found by substituting 



z 1 - Vl - Aw 

W = rrr Or Z = ; 

(1 + z) 2 1 + s/1 - Aw 

in f)5 .4[) . then setting 

sz 

W= (TT^ 

in (|5.6[) . The calculations used in obtaining the distribution and prob- 
ability generating function are the same as those used in obtaining the 
formula (I2.3[) in Griffiths 17f . There is a connection with a simple ran- 
dom walk on Z with transitions j — > j +1 with probability p and j — > j — 1 
with probability q< = 1 — p, when q > p. Let the number of steps to hit 
— 8, starting from 0, be £. Then £ has a probability generating function 
of 



2ps 

and i (£ + 8) has a probability generating function 



K(s) 



1 - VI - Apqs 
2p 
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A e (t) + 9 has the same distribution as the size-biased distribution of 
+ &)■, with probability generating function 

sK'(s) 



^A<>(t)( S ) — 



K'(l) 



identical to (|5.5p . In the random walk interpretation 9 is assumed to 
be an integer; however H(s) is infinitely divisible, so we use the same 
description for all 9 > 0. Another interpretation is that K(s) is the 
probability generating function of the total number of progeny in a 
Galton- Watson branching process with geometric offspring distribution 
qp k , k £ Z+, and extinction probability 1, beginning with 9 individu- 
als. See [11| Sections X.13 and XII. 5 for details of the random walk and 
branching process descriptions. An analogous calculation to (|5.3p which 



is included in Theorem 2.1 of 171 is that 



°(A e (s + t)=j\A< > (s)=i) 
i\ r(t + e)T(2j + 6) 



j)r(j + 9)r(i + 3+'9/ {1 z) 

x 2F1 {—i + 7 + 1, 2j + 9- 1 + j + 9; z), (5.7) 
where z — e . The jump rate from i — > j found from (|5.7[) is 

n ni±w|^) 2flH+ , +1 , 2 , +e;i+3 . +9;1) , 

\jj T(j + 6)T(i + j + 6) 

B(j + 9,i- j)- 1 [ l x 2 i +e - l (l - x)W-n- 2 dx 
3J Jo 

3) r(i-j)ro+e)r(2i+e-i) ~ ±5 z z ' • ■ ■ ' Q ^ 

2j+e (5.8) 



r(2j+8) / 1 



if i = 00. 



Bertoin [4| , [5| studies the genealogical structure of trees in an infinitely- 
many-alleles branching process model. In a limit from a large initial 
population size with rare mutations the genealogy is described by a 
continuous-state branching process in discrete time with an Inverse 
Gaussian reproduction law. We expect that there is a fascinating con- 
nection with the process {A e (t),t > 0}. A potential class of transition 
functions of Markov processes {qf.(t), t > 0} which are more general than 
subordinated processes and related to Bochner's characterization comes 
from replacing by p e n {t) by c n (t) described by (I3.7[) : however it is not 
clear that all such potential transition functions are positive, apart from 
those derived by subordination. 
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